###Load packages
library("rlang")
library("cli")
library("vctrs")
library("coda")
library("lattice")
library("MASS")
library("MCMCpack")
library("dplyr")
library("mcmcse")
library("polycor")
library("gridExtra")
library("ggmcmc")
library("stringr")
library("tibble")
library("readxl")
library("tidyverse")
library("dplyr")
library("ggplot2")

#Download reconnect_df.xlsx from Harvard Dataverse 

#Define the directory
setwd("")

###Cargar dataset
df <- read_excel("reconnect_df.xlsx")

###1. Preparing dataframe----

###Change column names
df <- df %>% rename(name = `Name of the Agency`,
                    country = Country,
                    pub_body = `Public Body`,
                    aut_recognition = `Autonomy recognition`,
                    term_office = `Term in office`,
                    term_renewal = `Term renewal`,
                    head_appoint = `Head appoint`,
                    president = President,
                    gov_board = `Governing Board`,
                    adv_board = `Advisory board`,
                    perf_survey = `Performing surveys`,
                    cons_decs = `Consulting decisions`,
                    pub_hear = `Public hearings`,
                    pub_enforce = `Publish enforcements`,
                    soc_media = `Disseminating info`,
                    pers_info = `Agency personnel info`,
                    cons_reps1 = `Consumer reps. in boards1`,
                    cons_reps2 = `Consumer reps. in boards2`,
                    cons_reps3 = `Consumer reps. in boards3`,
                    req_info = `Requesting information`,
                    pub_reports = `Publish reports`,
                    perf_research = `Research issues`, 
                    inspect_site = `Inspect on site`,
                    alert_sys = `Alert system`,
                    sanctions1 = `sanctioning capacity (4.22a)`,
                    sanctions2 = `sanc_number (4.22b)`,
                    sanctions3 = `sanc_business_closure (4.22c)`,
                    lit_cap = `Litigation capacity`,
                    norm_cap = `Normative capacity`,
                    formal_inv = `Formal investigations`,
                    conf_res1 = `Conflict resolution_capacity to promote`,
                    conf_res2 = `Conflict resolution_mechanism ADR in place`,
                    train_vulncons = `Vulnerable consumers_training 5.26c`,
                    educ_pub = `Educate public`,
                    online_cons = `Online consultation`,
                    online_recla = `Online reclamation`,
                    cit_den1 = `Citizens denounce companies1`,
                    cit_den2 = `Citizens denounce companies2`)

df$code <- as.factor(df$ID)
df$code <- dplyr::recode_factor(df$code, 
                                "AG1" = "DGCA-ESP", "AG2" = "CMA-GBR", 
                                "AG3" = "ACM-NLD", "AG4" = "NCA-NOR", 
                                "AG5" = "DGCCAFP-FRA", "AG6" = "SCA-SWE", 
                                "AG7" = "ICA-ITA", "AG8" = "FMSAHCCP-AUT",
                                "AG9" = "FMENCNSCP-DEU", "AG10" = "FCO-DEU",
                                "AG11" = "OCCP-POL", "AG12" = "CD-PRT",
                                "AG13" = "FPSESSE-BEL", "AG14" = "DCCA-DEN",
                                "AG15" = "DCO-DNK", "AG16" = "GDCP-GRC",
                                "AG17" = "CPTSA-EST", "AG18" = "NACP-ROU",
                                "AG19" = "FCCA-FIN", "AG20" = "CRPC-LVA",
                                "AG21" = "SCRPA-LTU", "AG22" = "MIT-CZE",
                                "AG23" = "CCPC-IRL", "AG24" = "CCPD-SVN",
                                "AG25" = "STI-SVK", "AG26" = "CPA-HUN",
                                "AG27" = "HCA-HUN", "AG28" = "CCP-BGR",
                                "AG29" = "SI-HRV", "AG30" = "CPS-CYP",
                                "AG31" = "CTIA-CZE", "AG32" = "MI-SVN",
                                "AG33" = "GDCM-ITA", "AG34" = "DCP-LUX",
                                "AG35" = "MCCAA-MLT", "AG36" = "CO-GRC",
                                "AG37" = "FCO-FIN", "AG38" = "DTIM-HRV")
table(df$code)

str(df$sanctions1)
str(df$sanctions2)
str(df$sanctions3)

df$sum_sanctions <- rowSums(df[, c("sanctions1", "sanctions2", "sanctions3")],
                            na.rm = TRUE)
table(df$sum_sanctions)

df$sum_sanctions <- ifelse(df$sum_sanctions == 2, 1, 
                           ifelse(df$sum_sanctions == 3, 2, 
                                  ifelse(df$sum_sanctions == 4, 3, 
                                         ifelse(df$sum_sanctions == 5, 4, df$sum_sanctions))))

###Descriptive variables
str(df$ID)
df$ID <- as.factor(df$ID)
table(df$ID)

table(df$country)
str(df$country)
df$country <- as.factor(df$country)

str(df$name)
df$name <- as.factor(df$name)

###Dimension 1: Autonomy
str(df$pub_body)
table(df$pub_body)
df$pub_body <- factor(df$pub_body, levels = c("0", "1", "2", "3"), 
                      ordered = TRUE)

str(df$norm_cap)
table(df$norm_cap)
df$norm_cap <- factor(df$norm_cap, levels = c("0", "1"), 
                      ordered = TRUE)

str(df$term_office)
table(df$term_office)

str(df$term_renewal)
table(df$term_renewal)
df$term_renewal <- factor(df$term_renewal, levels = c("0", "1", "2", "3"), 
                             ordered = TRUE)

str(df$head_appoint)
table(df$head_appoint)
df$head_appoint <- factor(df$head_appoint, levels = c("0", "1"), 
                          ordered = TRUE)

str(df$president)
table(df$president)
df$president <- factor(df$president, levels = c("0", "1"), 
                          ordered = TRUE)

str(df$gov_board)
table(df$gov_board)
df$gov_board <- factor(df$gov_board, levels = c("0", "1"), 
                       ordered = TRUE)

###Dimension 2: Accountability
str(df$adv_board)
table(df$adv_board)
df$adv_board <- factor(df$adv_board, levels = c("0", "1"), 
                       ordered = TRUE)

str(df$cons_reps3)
table(df$cons_reps3)
df$cons_reps3 <- factor(df$cons_reps3, levels = c("0", "1"), 
                        ordered = TRUE)

str(df$perf_survey)
table(df$perf_survey)
df$perf_survey <- factor(df$perf_survey, levels = c("0", "1"), 
                       ordered = TRUE)

str(df$cons_decs)
table(df$cons_decs)
df$cons_decs <- factor(df$cons_decs, levels = c("0", "1"), 
                         ordered = TRUE)

str(df$pub_hear)
table(df$pub_hear)
df$pub_hear <- factor(df$pub_hear, levels = c("0", "1"), 
                       ordered = TRUE)

str(df$req_info)
table(df$req_info)
df$req_info <- factor(df$req_info, levels = c("0", "1", "2"), 
                      ordered = TRUE)

str(df$pub_enforce)
table(df$pub_enforce)
df$pub_enforce <- factor(df$pub_enforce, levels = c("0", "1", "2"), 
                      ordered = TRUE)

str(df$soc_media)
table(df$soc_media)
df$soc_media <- factor(df$soc_media, levels = c("0", "1", "2", "3"), 
                         ordered = TRUE)

str(df$pers_info)
table(df$pers_info)
df$pers_info <- factor(df$pers_info, levels = c("0", "1", "2"), 
                       ordered = TRUE)

str(df$pub_reports)
table(df$pub_reports)
df$pub_reports <- factor(df$pub_reports, levels = c("0", "1"), 
                      ordered = TRUE)

###Dimension 3: Authorities' enforcement
str(df$perf_research)
table(df$perf_research)
df$perf_research <- factor(df$perf_research, levels = c("0", "1"), 
                         ordered = TRUE)

str(df$formal_inv)
table(df$formal_inv)
df$formal_inv <- factor(df$formal_inv, levels = c("0", "1"), 
                        ordered = TRUE)

str(df$inspect_site)
table(df$inspect_site)
df$inspect_site <- factor(df$inspect_site, levels = c("0", "1"), 
                           ordered = TRUE)

str(df$alert_sys)
table(df$alert_sys)
df$alert_sys <- factor(df$alert_sys, levels = c("0", "1"), 
                          ordered = TRUE)

str(df$sum_sanctions)
table(df$sum_sanctions)

str(df$lit_cap)
table(df$lit_cap)
df$lit_cap <- factor(df$lit_cap, levels = c("0", "1"), 
                       ordered = TRUE)

str(df$conf_res1)
table(df$conf_res1)
df$conf_res1 <- factor(df$conf_res1, levels = c("0", "1"), 
                        ordered = TRUE)

str(df$conf_res2)
table(df$conf_res2)
df$conf_res2 <- factor(df$conf_res2, levels = c("0", "1", "2", "3"), 
                       ordered = TRUE)

###Dimension 4: consumers' awareness
str(df$train_vulncons)
table(df$train_vulncons)
df$train_vulncons <- factor(df$train_vulncons, levels = c("0", "1"), 
                       ordered = TRUE)

str(df$educ_pub)
table(df$educ_pub)
df$educ_pub <- factor(df$educ_pub, levels = c("0", "1"), 
                            ordered = TRUE)

str(df$online_cons)
table(df$online_cons)
df$online_cons <- factor(df$online_cons, levels = c("0", "1"), 
                      ordered = TRUE)

str(df$online_recla)
table(df$online_recla)
df$online_recla <- factor(df$online_recla, levels = c("0", "1"), 
                         ordered = TRUE)

str(df$cit_den1)
table(df$cit_den1)
df$cit_den1 <- factor(df$cit_den1, levels = c("0", "1"), 
                          ordered = TRUE)

str(df$cit_den2)
table(df$cit_den2)
df$cit_den2 <- factor(df$cit_den2, levels = c("0", "1"), 
                      ordered = TRUE)

##Plot agencies by institutional design (Figure 1)
table(df$pub_body)
df$pub_body_cat <- dplyr::recode_factor(df$pub_body, 
                                        "0" = "Unit within Ministry",
                                        "1" = "Ministry",
                                        "2" = "Unit within Agency",
                                        "3" = "Public Agency"
                                        ) 

table(df$pub_body_cat)

df_counts <- df %>%
  count(pub_body_cat)

ggplot(df_counts, aes(x = pub_body_cat, y = n)) +
  geom_bar(stat = "identity", fill = "gray60", width = 0.6) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(
    x = NULL,
    y = "Observations",
    ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  scale_y_continuous(
    breaks = seq(0, 14, by = 2),
    limits = c(0, 14),
    expand = c(0, 0)
  )

ggsave("fig1.png",
       width = 9, height = 6, dpi = 300)

####2. Bayesian factorial analysis----
###Running the chain----
###Paso 2. Priors 
#### The priors serve as a reference point for estimating the Bayesian factors.
# They are calculated based on the mean of each item used in the Bayesian factors,
# drawing on data from two randomly selected agencies.

## The selected agencies, along with their values across the dimension items,
## are documented in the "sample_agencies" file in the repository.

library("openxlsx")

attach(df)

###Priors: Autonomy dimension
D1.prior <- c(0,0,0,0,0,0,1,
              0,0,0,0,0,0,1) 
dim(D1.prior) <- c(7,2)

###Priors: Accountability dimension
D2.prior <- c(1,1,1,1,1,1,2,3,0,1,
              1,1,1,1,1,1,0,1,1,1)
dim(D2.prior) <- c(10,2)

###Priors: Enforcement dimension
D3.prior <- c(1,1,1,3,1,1,0,3,
              1,0,0,0,0,0,0,3)
dim(D3.prior) <- c(8,2)

###Priors: Consumers' Awareness
D4.prior <- c(1,1,1,1,1,1,
              0,1,1,0,0,0)
dim(D4.prior) <- c(6,2)

###Paso 3. Chain
###Autonomy----
D1.mcmc <- MCMCmixfactanal(~pub_body+
                             term_office+
                             term_renewal+
                             head_appoint+
                             president+
                             gov_board+
                             norm_cap,
                           lambda.constraints = 
                             list(),
                           seed = 1090,
                           factors = 1,
                           data= df,
                           burnin = 10000,
                           mcmc = 1000000,
                           thin=100, 
                           verbose = 10000,
                           L0 = 0.25,
                           store.lambda = TRUE,
                           store.scores = TRUE,
                           tune = .25, 
                           l0 = D1.prior)

save.image() # Save chain

view(D1.mcmc)


## Preparing for analysis: extracting labels for later 
# Dimension names
dimnames.1 <- attr(D1.mcmc, which = "dimnames")
names.1 <- dimnames.1[[2]]
view(names.1)

# Filter
lambda1 <-  grep("^Lambda", names.1, value = TRUE) # All lambdas
view(lambda1)

lambda1.1 <- grep("1$", lambda1, value = TRUE)  # Lambda 1 (negative item difficulty parameter)
view(lambda1.1)

lambda1.2 <- grep("2$", lambda1, value = TRUE) # Lambda 2 (factor loadings or the item discrimination parameters)
view(lambda1.2)

phi1 <-  grep("phi", names.1, value = TRUE) # Phi
view(phi1)

# Create labels
lambda1.1.label <- lambda1.1 %>% str_sub(start = 7, end = -3)%>% paste0(" (lambda 1)")
view(lambda1.1.label)

lambda1.2.label <- lambda1.2 %>% str_sub(start = 7, end = -3) %>% paste0(" (lambda 2)")
view(lambda1.2.label)

agencies.label1  <- phi1 %>% str_sub(start = 5, end = -3) %>% paste0(" (agencies)")
view(agencies.label1)

agencies.label.clean1  <- phi1 %>% str_sub(start = 5, end = -3)
view(agencies.label.clean1)

# Create dataframe
lambda1.1.df <- data.frame(
  Parameter = lambda1.1,
  Label = lambda1.1.label)
view(lambda1.1.df)

lambda1.2.df <- data.frame(
  Parameter = lambda1.2,
  Label = lambda1.2.label)
view(lambda1.2.df)

agencies1.df <- data.frame(
  Parameter = phi1,
  Label = agencies.label1)
view(agencies1.df)

agencies.clean1.df <- data.frame(
  Parameter = phi1,
  Label = agencies.label.clean1)
view(agencies.clean1.df)

total1.df <- rbind(agencies1.df, agencies.clean1.df, lambda1.1.df, 
                   lambda1.2.df) 
view(total1.df)

# Convert mcmc object
sample.post.converted1 <- ggs(as.mcmc.list(D1.mcmc), par_labels = total1.df)
view(sample.post.converted1)

sample.post.agencies1 <- ggs(as.mcmc.list(D1.mcmc), par_labels = agencies.clean1.df, family = "^phi")
view(sample.post.agencies1)

sample.post.items1 <- ggs(as.mcmc.list(D1.mcmc), par_labels = agencies.clean1.df, family = "Lambda")
view(sample.post.items1)

sample.post.items1 <- sample.post.items1 %>%
  filter(str_detect(Parameter, "\\.1") | Parameter == "Lambdaterm_office.2") %>%
  filter(str_detect(Parameter, "Lambda")) %>%
  mutate(Parameter = str_remove(Parameter, "\\.1"),
         Parameter = str_remove(Parameter, "Lambda"))

table(sample.post.items1$Parameter)

## Key information
# Summary
summary <- summary(D1.mcmc)
scores1.df <- summary[["statistics"]] %>%
  as.data.frame() 

scores1 <- scores1.df %>%
  mutate(Agency = row.names(scores1.df)) %>%
  filter(str_detect(Agency, "^phi"))  %>%
  mutate(Agency = str_remove(Agency, "phi.")) %>%
  mutate(Agency = substr(Agency, 1, nchar(Agency) - 2))

#Change direction of factorial weights (Means)
scores1$Mean <- -scores1$Mean

scores.items1 <- scores1.df %>%
  mutate(Item = row.names(scores1.df)) %>%
  filter(str_detect(Item, ".1")) %>%
  filter(str_detect(Item, "Lambda")) %>%
  mutate(Item = str_remove(Item, ".1"),
         Item = str_remove(Item, "Lambda"))

##Change the direction of factorial weights
scores.items1$Mean <- -scores.items1$Mean

###Plot scores by agencies

str(sample.post.agencies1$value)
summary(sample.post.agencies1$value)
sample.post.agencies1$value <- -sample.post.agencies1$value

table(sample.post.agencies1$ParameterOriginal)
sample.post.agencies1$ID <- as.character(sample.post.agencies1$ParameterOriginal)
table(sample.post.agencies1$ID)
sample.post.agencies1$ID <- gsub("phi\\.([0-9]+)\\.2", "AG\\1", 
                                 sample.post.agencies1$ID)
sample.post.agencies1$ID <- as.factor(sample.post.agencies1$ID)

sample.post.agencies1$ID <- factor(sample.post.agencies1$ID, 
                                   levels = paste0("AG", 1:38))

sample.post.agencies1 <- merge(sample.post.agencies1, df[, c("ID", "code")], 
                               by = "ID", all.x = TRUE)
str(sample.post.agencies1$code)
str(sample.post.agencies1$Parameter)

sample.post.agencies1$Parameter <- sample.post.agencies1$code

# Plot items score (Figure 2)
##Change the direction of factorial weights
sample.post.items1$value <- -sample.post.items1$value

ggs_caterpillar(sample.post.items1, line = 0) +
  scale_x_continuous(name =  "Autonomy",
                     breaks = seq(-6, 7, 1),  
                     limits = c(-6, 7)) +
  scale_y_discrete(name = "Items",
                   labels=c("pub_body" = "Type of public\nbody",
                            "president" = "Non-executive\npresident",
                            "term_office.2" = "Agency head\nterm of office",
                            "term_renewal" = "Agency head\nextension",
                            "norm_cap" = "Normative\ncapacity",
                            "gov_board" = "Government\nboard",
                            "head_appoint" = "Agency head\nappointment")) + 
  labs(y = "") + theme_linedraw() +
  theme(axis.title.x  = element_text(size = 15, face = "italic"),
        axis.title.y  = element_text(size = 14),           
        axis.text.y   = element_text(size = 13),           
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("D1items.png",
       width = 9, height = 6, dpi = 300)

###Paso 3. Chain
###Accountability----
D2.mcmc <- MCMCordfactanal(~adv_board+
                             cons_reps3+
                             perf_survey+
                             cons_decs+
                             pub_hear+
                             req_info+
                             pub_enforce+
                             soc_media+
                             pers_info+
                             pub_reports,
                           lambda.constraints = 
                             list(),
                           seed = 1090,
                           factors = 1,
                           data= df,
                           burnin = 10000,
                           mcmc = 1000000,
                           thin=100, 
                           verbose = 10000,
                           L0 = 0.25,
                           store.lambda = TRUE,
                           store.scores = TRUE,
                           tune = .25, 
                           l0 = D2.prior)

save.image() # Save chain

view(D2.mcmc)

## Preparing for analysis: extracting labels for later 

# Dimension names
dimnames.2 <- attr(D2.mcmc, which = "dimnames")
names.2 <- dimnames.2[[2]]

# Filter
lambda2 <-  grep("^Lambda", names.2, value = TRUE) # All lambdas
view(lambda2)

lambda2.1 <- grep("1$", lambda2, value = TRUE)  # Lambda 1 (negative item difficulty parameter)
view(lambda2.1)

lambda2.2 <- grep("2$", lambda2, value = TRUE) # Lambda 2 (factor loadings or the item discrimination parameters)
view(lambda2.2)

phi2 <-  grep("phi", names.2, value = TRUE) # Phi
view(phi2)

# Create labels
lambda2.1.label <- lambda2.1 %>% str_sub(start = 7, end = -3)%>% paste0(" (lambda 1)")
view(lambda2.1.label)

lambda2.2.label <- lambda2.2 %>% str_sub(start = 7, end = -3) %>% paste0(" (lambda 2)")
view(lambda2.2.label)

agencies.label2  <- phi2 %>% str_sub(start = 5, end = -3) %>% paste0(" (agencies)")
view(agencies.label2)

agencies.label.clean2  <- phi2 %>% str_sub(start = 5, end = -3)
view(agencies.label.clean2)

# Create dataframe
lambda2.1.df <- data.frame(
  Parameter = lambda2.1,
  Label = lambda2.1.label)
view(lambda2.1.df)

lambda2.2.df <- data.frame(
  Parameter = lambda2.2,
  Label = lambda2.2.label)
view(lambda2.2.df)

agencies2.df <- data.frame(
  Parameter = phi2,
  Label = agencies.label2)
view(agencies2.df)

agencies.clean2.df <- data.frame(
  Parameter = phi2,
  Label = agencies.label.clean2)
view(agencies.clean2.df)

total2.df <- rbind(agencies2.df, agencies.clean2.df, lambda2.1.df, 
                   lambda2.2.df)
view(total2.df)

# Convert mcmc object
sample.post.converted2 <- ggs(as.mcmc.list(D2.mcmc), par_labels = total2.df)
view(sample.post.converted2)

sample.post.agencies2 <- ggs(as.mcmc.list(D2.mcmc), par_labels = agencies.clean2.df, family = "^phi")
view(sample.post.agencies2)

sample.post.items2 <- ggs(as.mcmc.list(D2.mcmc), par_labels = agencies.clean2.df, family = "Lambda")
view(sample.post.items2)

sample.post.items2 <- sample.post.items2%>%
  filter(str_detect(Parameter, ".2")) %>%
  filter(str_detect(Parameter, "Lambda")) %>%
  mutate(Parameter = str_remove(Parameter, ".2"),
         Parameter = str_remove(Parameter, ".2"), 
         Parameter = str_remove(Parameter, "Lambda"))

table(sample.post.items2$Parameter)

## Key information
# Summary
summary <- summary(D2.mcmc)

scores2.df <- summary[["statistics"]] %>%
  as.data.frame() 

scores2 <- scores2.df %>%
  mutate(Agency = row.names(scores2.df)) %>%
  filter(str_detect(Agency, "^phi"))  %>%
  mutate(Agency = str_remove(Agency, "phi.")) %>%
  mutate(Agency = substr(Agency, 1, nchar(Agency) - 2))

scores.items2 <- scores2.df %>%
  mutate(Item = row.names(scores2.df)) %>%
  filter(str_detect(Item, ".2")) %>%
  filter(str_detect(Item, "Lambda")) %>%
  mutate(Item = str_remove(Item, ".2"),
         Item = str_remove(Item, "Lambda"))

table(sample.post.agencies2$ParameterOriginal)
sample.post.agencies2$ID <- as.character(sample.post.agencies2$ParameterOriginal)
table(sample.post.agencies2$ID)
sample.post.agencies2$ID <- gsub("phi\\.([0-9]+)\\.2", "AG\\1", 
                                 sample.post.agencies2$ID)
sample.post.agencies2$ID <- as.factor(sample.post.agencies2$ID)

sample.post.agencies2$ID <- factor(sample.post.agencies2$ID, 
                                   levels = paste0("AG", 1:38))

sample.post.agencies2 <- merge(sample.post.agencies2, df[, c("ID", "code")], 
                               by = "ID", all.x = TRUE)
str(sample.post.agencies2$code)
str(sample.post.agencies2$Parameter)

sample.post.agencies2$Parameter <- sample.post.agencies2$code

###Plot (Figure 3)

ggs_caterpillar(sample.post.items2, line = 0) +
  scale_x_continuous(name =  "Accountability",
                     breaks = seq(-3, 7, 1),  
                     limits = c(-3, 7)) +
  scale_y_discrete(name = "Items",
                   labels=c("adv_board" = "Advisory\nboard",
                            "perf_survey" = "Performing\nsurveys",
                            "cons_decs" = "Consulting\ndecisions",
                            "pub_hear" = "Public\nhearings",
                            "pub_enforce" = "Publishing resolutions\non the website",
                            "soc_media" = "Use of\nsocial media",
                            "pers_info" = "Agency personnel\ninfo.",
                            "cons_reps3" = "Consumer reps.\non boards",
                            "req_info" = "Request info.\non the website",
                            "pub_reports" = "Reports\npublished")) + 
  labs(y = "") + theme_linedraw() +
  theme(axis.title.x  = element_text(size = 15, face = "italic"),
        axis.title.y  = element_text(size = 14),           
        axis.text.y   = element_text(size = 13),           
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("D2items.png", width = 9, height = 6, dpi = 300)

###Paso 3. Chain
###Authorities' Enforcement----
D3.mcmc <- MCMCmixfactanal(~perf_research+
                             inspect_site+
                             alert_sys+
                             sum_sanctions+
                             lit_cap+
                             formal_inv+
                             conf_res1+
                             conf_res2,
                           lambda.constraints = 
                             list(),
                           seed = 1090,
                           factors = 1,
                           data= df,
                           burnin = 10000,
                           mcmc = 1000000,
                           thin=100, 
                           verbose = 10000,
                           L0 = 0.25,
                           store.lambda = TRUE,
                           store.scores = TRUE,
                           tune = .25, 
                           l0 = D3.prior)

save.image() # Save chain

view(D3.mcmc)

## Preparing for analysis: extracting labels for later 

# Dimension names
dimnames.3 <- attr(D3.mcmc, which = "dimnames")
names.3 <- dimnames.3[[2]]

# Filter
lambda3 <-  grep("^Lambda", names.3, value = TRUE) # All lambdas
view(lambda3)

lambda3.1 <- grep("1$", lambda3, value = TRUE)  # Lambda 1 (negative item difficulty parameter)
view(lambda3.1)

lambda3.2 <- grep("2$", lambda3, value = TRUE) # Lambda 2 (factor loadings or the item discrimination parameters)
view(lambda3.2)

phi3 <-  grep("phi", names.3, value = TRUE) # Phi
view(phi3)

# Create labels
lambda3.1.label <- lambda3.1 %>% str_sub(start = 7, end = -3)%>% paste0(" (lambda 1)")
view(lambda3.1.label)

lambda3.2.label <- lambda3.2 %>% str_sub(start = 7, end = -3) %>% paste0(" (lambda 2)")
view(lambda3.2.label)

agencies.label3 <- phi3 %>% str_sub(start = 5, end = -3) %>% paste0(" (agencies)")
view(agencies.label3)

agencies.label.clean3 <- phi3 %>% str_sub(start = 5, end = -3)
view(agencies.label.clean3)

# Create dataframe
lambda3.1.df <- data.frame(
  Parameter = lambda3.1,
  Label = lambda3.1.label)

lambda3.2.df <- data.frame(
  Parameter = lambda3.2,
  Label = lambda3.2.label)

agencies3.df <- data.frame(
  Parameter = phi3,
  Label = agencies.label3)

agencies.clean3.df <- data.frame(Parameter = phi3,
                                 Label = agencies.label.clean3)

total3.df <- rbind(agencies3.df, agencies.clean3.df, lambda3.1.df, 
                   lambda3.2.df)#, gamma3.df)
view(total3.df)

# Convert mcmc object
sample.post.converted3 <- ggs(as.mcmc.list(D3.mcmc), par_labels = total3.df)
sample.post.agencies3 <- ggs(as.mcmc.list(D3.mcmc), par_labels = agencies.clean3.df, family = "^phi")
sample.post.items3 <- ggs(as.mcmc.list(D3.mcmc), par_labels = agencies.clean3.df, family = "Lambda")
table(sample.post.items3$Parameter)

sample.post.items3 <- sample.post.items3 %>%
  filter(str_detect(Parameter, "\\.2")) %>% 
  filter(str_detect(Parameter, "Lambda")) %>%  
  mutate(Parameter = str_remove(Parameter, "Lambda"),  
         Parameter = str_remove(Parameter, "\\.2")) 

## Key information

# Summary
summary <- summary(D3.mcmc)

scores3.df <- summary[["statistics"]] %>%
  as.data.frame() 

scores3 <- scores3.df %>%
  mutate(Agency = row.names(scores3.df)) %>%
  filter(str_detect(Agency, "^phi"))  %>%
  mutate(Agency = str_remove(Agency, "phi.")) %>%
  mutate(Agency = substr(Agency, 1, nchar(Agency) - 2))

#Change direction of factorial weights (Means)
scores3$Mean <- -scores3$Mean

scores.items3 <- scores3.df %>%
  mutate(Item = row.names(scores3.df)) %>%
  filter(str_detect(Item, ".2")) %>%
  filter(str_detect(Item, "Lambda")) %>%
  mutate(Item = str_remove(Item, ".2"),
         Item = str_remove(Item, "Lambda"))

##Change the direction of factorial weights
scores.items3$Mean <- -scores.items3$Mean

table(sample.post.agencies3$ParameterOriginal)
sample.post.agencies3$ID <- as.character(sample.post.agencies3$ParameterOriginal)
table(sample.post.agencies3$ID)
sample.post.agencies3$ID <- gsub("phi\\.([0-9]+)\\.2", "AG\\1", 
                                 sample.post.agencies3$ID)
sample.post.agencies3$ID <- as.factor(sample.post.agencies3$ID)

sample.post.agencies3$ID <- factor(sample.post.agencies3$ID, 
                                   levels = paste0("AG", 1:38))

sample.post.agencies3 <- merge(sample.post.agencies3, df[, c("ID", "code")], 
                               by = "ID", all.x = TRUE)
str(sample.post.agencies3$code)
str(sample.post.agencies3$Parameter)

sample.post.agencies3$Parameter <- sample.post.agencies3$code

###Plot (Figure 4)

sample.post.items3$value <- -sample.post.items3$value

ggs_caterpillar(sample.post.items3, line = 0) +
  scale_x_continuous(name =  "Authorities' Enforcement",
                     breaks = seq(-3, 7, 1),  
                     limits = c(-5, 7)) +
  scale_y_discrete(name = "Items",
                   labels=c("perf_research" = "Perfoming\nresearch",
                            "inspect_site" = "Inspections\non site",
                            "alert_sys" = "Alert system",
                            "sum_sanctions" = "Summ. sanction\nindex",
                            "lit_cap" = "Litigation\ncapacity",
                            "formal_inv" = "Formal\ninvestigations",
                            "conf_res1" = "Conf. res.:\ncapacity to promote",
                            "conf_res2" = "Conf. res.:ADR\nmechanism in place")) + 
  labs(y = "") + theme_linedraw() +
  theme(axis.title.x  = element_text(size = 15, face = "italic"),
        axis.title.y  = element_text(size = 14),           
        axis.text.y   = element_text(size = 13),           
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("D3items.png", width = 9, height = 6, dpi = 300)

###Consumers' awareness----
D4.mcmc <- MCMCordfactanal(~train_vulncons+
                             educ_pub+
                             online_cons+
                             online_recla+
                             cit_den1+
                             cit_den2,
                           lambda.constraints = 
                             list(),
                           seed = 1090,
                           factors = 1,
                           data= df,
                           burnin = 10000,
                           mcmc = 1000000,
                           thin=100, 
                           verbose = 10000,
                           L0 = 0.25,
                           store.lambda = TRUE,
                           store.scores = TRUE,
                           tune = .25, 
                           l0 = D4.prior)

save.image() # Save chain

view(D4.mcmc)

## Preparing for analysis: extracting labels for later 
# Dimension names
dimnames.4 <- attr(D4.mcmc, which = "dimnames")
names.4 <- dimnames.4[[2]]
view(names.4)

# Filter
lambda4 <-  grep("^Lambda", names.4, value = TRUE) # All lambdas
view(lambda4)

lambda4.1 <- grep("1$", lambda4, value = TRUE)  # Lambda 1 (negative item difficulty parameter)
view(lambda4.1)

lambda4.2 <- grep("2$", lambda4, value = TRUE) # Lambda 2 (factor loadings or the item discrimination parameters)
view(lambda4.2)

phi4 <-  grep("phi", names.4, value = TRUE) # Phi
view(phi4)

# Create labels
lambda4.1.label <- lambda4.1 %>% str_sub(start = 7, end = -3)%>% paste0(" (lambda 1)")
view(lambda4.1.label)

lambda4.2.label <- lambda4.2 %>% str_sub(start = 7, end = -3) %>% paste0(" (lambda 2)")
view(lambda4.2.label)

agencies.label4  <- phi4 %>% str_sub(start = 5, end = -3) %>% paste0(" (agencies)")
view(agencies.label4)

agencies.label.clean4  <- phi4 %>% str_sub(start = 5, end = -3)
view(agencies.label.clean4)

# Create dataframe
lambda4.1.df <- data.frame(
  Parameter = lambda4.1,
  Label = lambda4.1.label)
view(lambda4.1.df)

lambda4.2.df <- data.frame(Parameter = lambda4.2,
                           Label = lambda4.2.label)
view(lambda4.2.df)

agencies4.df <- data.frame(
  Parameter = phi4,
  Label = agencies.label4)
view(agencies4.df)

agencies.clean4.df <- data.frame(
  Parameter = phi4,
  Label = agencies.label.clean4)
view(agencies.clean4.df)

total4.df <- rbind(agencies4.df, agencies.clean4.df, lambda4.1.df, 
                   lambda4.2.df)
view(total4.df)

# Convert mcmc object
sample.post.converted4 <- ggs(as.mcmc.list(D4.mcmc), par_labels = total4.df)
view(sample.post.converted4)

sample.post.agencies4 <- ggs(as.mcmc.list(D4.mcmc), par_labels = agencies.clean4.df, family = "^phi")
view(sample.post.agencies4)

sample.post.items4 <- ggs(as.mcmc.list(D4.mcmc), par_labels = agencies.clean4.df, family = "Lambda")
view(sample.post.items4)
table(sample.post.items4$Parameter)

sample.post.items4 <- sample.post.items4%>%
  filter(str_detect(Parameter, "\\.2")) %>%
  filter(str_detect(Parameter, "Lambda")) %>%
  mutate(Parameter = str_remove(Parameter, "Lambda"),
         Parameter = str_remove(Parameter, "\\.2"))

## Key information
# Summary
summary <- summary(D4.mcmc)

scores4.df <- summary[["statistics"]] %>%
  as.data.frame() 

scores4 <- scores4.df %>%
  mutate(Agency = row.names(scores4.df)) %>%
  filter(str_detect(Agency, "^phi"))  %>%
  mutate(Agency = str_remove(Agency, "phi.")) %>%
  mutate(Agency = substr(Agency, 1, nchar(Agency) - 2))

scores.items4 <- scores4.df %>%
  mutate(Item = row.names(scores4.df)) %>%
  filter(str_detect(Item, ".2")) %>%
  filter(str_detect(Item, "Lambda")) %>%
  mutate(Item = str_remove(Item, ".2"),
         Item = str_remove(Item, "Lambda"))

table(sample.post.agencies4$ParameterOriginal)
sample.post.agencies4$ID <- as.character(sample.post.agencies4$ParameterOriginal)
table(sample.post.agencies4$ID)
sample.post.agencies4$ID <- gsub("phi\\.([0-9]+)\\.2", "AG\\1", 
                                 sample.post.agencies4$ID)
sample.post.agencies4$ID <- as.factor(sample.post.agencies4$ID)

sample.post.agencies4$ID <- factor(sample.post.agencies4$ID, 
                                   levels = paste0("AG", 1:38))

sample.post.agencies4 <- merge(sample.post.agencies4, df[, c("ID", "code")], 
                               by = "ID", all.x = TRUE)
str(sample.post.agencies4$code)
str(sample.post.agencies4$Parameter)

sample.post.agencies4$Parameter <- sample.post.agencies4$code

##Plot (Figure 5)

ggs_caterpillar(sample.post.items4, line = 0) +
  scale_x_continuous(name =  "Consumers awareness",
                     breaks = seq(-3, 7, 1),  
                     limits = c(-3, 7)) +
  scale_y_discrete(name = "Items",
                   labels=c("train_vulncons" = "Training offer about\nvulnerable consumers",
                            "educ_pub" = "Programs to\neducate consumers",
                            "online_cons" = "Online consultations\navailable",
                            "online_recla" = "Online complaints",
                            "cit_den1" = "Online reporting\nof firms",
                            "cit_den2" = "Online grievance-\ninformal complaint")) + 
  labs(y = "") + theme_linedraw() +
  theme(axis.title.x  = element_text(size = 15, face = "italic"),
        axis.title.y  = element_text(size = 14),           
        axis.text.y   = element_text(size = 13),           
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("D4items.png", width = 9, height = 5, dpi = 300)

###3. Scatter plots----
###Ceate dataframe with agencies' scores
df$ID2 <- 1:38
str(df$ID2)
df$ID2 <- as.factor(df$ID2)
table(df$ID2)

scores1$ID2 <- 1:38
str(scores1$ID2)
scores1$ID2 <- as.factor(scores1$ID2)
table(scores1$ID2)

scores2$ID2 <- 1:38
str(scores2$ID2)
scores2$ID2 <- as.factor(scores2$ID2)
table(scores2$ID2)

scores3$ID2 <- 1:38
str(scores3$ID2)
scores3$ID2 <- as.factor(scores3$ID2)
table(scores3$ID2)

scores4$ID2 <- 1:38
str(scores4$ID2)
scores4$ID2 <- as.factor(scores4$ID2)
table(scores4$ID2)

dfgen_scatter <- data.frame(df$ID2,df$name,df$country,df$code)
dfgen_scatter <- rename(dfgen_scatter, ID2 = df.ID2)
dfgen_scatter <- rename(dfgen_scatter, name = df.name)
dfgen_scatter <- rename(dfgen_scatter, country = df.country)
dfgen_scatter <- rename(dfgen_scatter, code = df.code)

scores1 <- rename(scores1, Mean.D1 = Mean)

dfgen_scatter <- merge(dfgen_scatter, scores1[, c("ID2", "Mean.D1")], 
                    by = "ID2")

scores2 <- rename(scores2, Mean.D2 = Mean)

dfgen_scatter <- merge(dfgen_scatter, scores2[, c("ID2", "Mean.D2")], 
                    by = "ID2")

scores3 <- rename(scores3, Mean.D3 = Mean)

dfgen_scatter <- merge(dfgen_scatter, scores3[, c("ID2", "Mean.D3")], 
                    by = "ID2")

scores4 <- rename(scores4, Mean.D4 = Mean)

dfgen_scatter <- merge(dfgen_scatter, scores4[, c("ID2", "Mean.D4")], 
                       by = "ID2")

table(dfgen_scatter$ID2)
levels(dfgen_scatter$ID2)

dfgen_scatter <- dfgen_scatter %>%
  arrange(ID2)

###Scatter plots
library("patchwork")
library("RColorBrewer")
library("viridisLite")
library("viridis")
library("ggrepel")

##Accountability by Agency Autonomy (Figure 6)
ggplot(dfgen_scatter, aes(x = Mean.D1, y = Mean.D2)) +
  geom_point(size = 2) +
  geom_text_repel(aes(label = code), 
                  box.padding = 0.175, 
                  point.padding = 0.1,
                  segment.color = "transparent",
                  show.legend = FALSE) +
  xlab("Autonomy") +
  ylab("Accountability") +
  xlim(-2.5, 2.5) +
  ylim(-2.5, 2) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "#FF9999", size = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "#FF9999", size = 1) +
  theme(
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    axis.text.x  = element_text(size = 13),
    axis.text.y  = element_text(size = 13)
  )

ggsave("accounxaut.png", width = 13, height = 8, dpi = 300)

###Consumers Awareness by Enforcement (Figure 7)
ggplot(dfgen_scatter, aes(x = Mean.D3, y = Mean.D4)) +
  geom_point(size = 2) +
  geom_text_repel(aes(label = code), 
                  box.padding = 0.125, 
                  point.padding = 0.5,
                  segment.color = "transparent",
                  show.legend = FALSE) +
  xlab("Auth. enforcement") +
  ylab("Consumers awareness") +
  xlim(-3, 2) +
  ylim(-2, 2) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "#FF9999", size = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "#FF9999", size = 1) +
  theme(
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    axis.text.x  = element_text(size = 13),
    axis.text.y  = element_text(size = 13)
  )

ggsave("consxenforce.png", width = 13, height = 8, dpi = 300)

####4. Cluster analysis----
library("dendextend")

df_dend <- dfgen_scatter
df_dend$code_2 <- as.character(df_dend$code)
str(df_dend$code_2)
colnames(df_dend)

df.dist <- dist(df_dend %>% dplyr::select(-c(ID2, name, country, code, code_2)),
                method = "euclidean")
df.hclust <- hclust(df.dist, method = "ward.D")

ordered_labels <- df_dend$code_2[df.hclust$order]  
df.dend <- as.dendrogram(df.hclust) %>%
  set("labels", as.character(ordered_labels))  

#Dendogram plot
plot(df.dend, cex = 0.4)

###Constructing cluster membership variable based on dendogram results

df_dend$dend.clust[df_dend$code_2=="MIT-CZE"|df_dend$code_2=="CCPD-SVN"|
                   df_dend$code_2=="GDCM-ITA"|df_dend$code_2=="DTIM-HRV"|
                   df_dend$code_2=="FCO-DEU"|df_dend$code_2=="DCCA-DEN"|
                   df_dend$code_2=="HCA-HUN"|df_dend$code_2=="FMSAHCCP-AUT"|
                   df_dend$code_2=="FMENCNSCP-DEU"|
                   df_dend$code_2=="DCP-LUX"] <- "Cluster 1"

df_dend$dend.clust[df_dend$code_2=="SCA-SWE"|df_dend$code_2=="CD-PRT"|
                         df_dend$code_2=="CRPC-LVA"|df_dend$code_2=="DGCA-ESP"|
                         df_dend$code_2=="CPS-CYP"|df_dend$code_2=="GDCP-GRC"|
                         df_dend$code_2=="DGCCAFP-FRA"|
                         df_dend$code_2=="FPSESSE-BEL"] <- "Cluster 2"

df_dend$dend.clust[df_dend$code_2=="STI-SVK"|df_dend$code_2=="MI-SVN"|
                   df_dend$code_2=="CPTSA-EST"|df_dend$code_2=="CPA-HUN"|
                   df_dend$code_2=="CTIA-CZE"] <- "Cluster 3"

df_dend$dend.clust[df_dend$code_2=="CO-GRC"|df_dend$code_2=="NACP-ROU"|
                     df_dend$code_2=="FCCA-FIN"|df_dend$code_2=="SI-HRV"|
                     df_dend$code_2=="FCO-FIN"|df_dend$code_2=="SCRPA-LTU"|
                     df_dend$code_2=="CCP-BGR"|df_dend$code_2=="MCCAA-MLT"|
                     df_dend$code_2=="OCCP-POL"|df_dend$code_2=="ACM-NLD"|
                     df_dend$code_2=="CCPC-IRL"|df_dend$code_2=="ICA-ITA"|
                     df_dend$code_2=="NCA-NOR"|df_dend$code_2=="CMA-GBR"|
                     df_dend$code_2=="DCO-DNK"] <- "Cluster 4"

table(df_dend$dend.clust)

##Plot the Ward.D cluters (Figure 8)
clust_mean <- df_dend%>%
  group_by(dend.clust) %>%
  summarise(
    Autonomy = mean(Mean.D1, na.rm = TRUE),  
    Accountability = mean(Mean.D2, na.rm = TRUE),  
    Enforcement = mean(Mean.D3, na.rm = TRUE),  
    C.Awareness = mean(Mean.D4, na.rm = TRUE),
    .groups = "drop"
  )

clust_long <- gather(clust_mean, key = "Variable", 
                      value = "Value", -dend.clust)

table(clust_long$Variable)
clust_long$Variable <- dplyr::recode_factor(clust_long$Variable, 
                                            "Autonomy" = "Autonomy",
                                            "Accountability" = "Accountability",
                                            "Enforcement" = "Enforcement",
                                            "C.Awareness" = "Consumers\nAwareness")

ggplot(clust_long, aes(x = Variable, y = Value, fill = "Group")) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.7, color = "#FFB6C1") +
  labs(x = "Dimensions", y = "Means of clusters' dimensions") +
  facet_wrap(~dend.clust, ncol = 1) +
  theme_minimal() + 
  ylim(c(-1.5, 1.5)) +
  theme(
    legend.position = "none",
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    axis.text.x  = element_text(size = 13),
    axis.text.y  = element_text(size = 13),
    strip.text = element_text(size = 14)  
  )

ggsave("clusters_wardD_prior3.png",
       width = 13, height = 8, dpi = 300)


###5. Regression analysis----
###Create binary variable for each cluster
table(df_dend$dend.clust)
str(df_dend$dend.clust)

df_dend$cluster1 <- dplyr::recode_factor(df_dend$dend.clust,
                                          "Cluster 1" = "1",
                                          .default = "0")
table(df_dend$cluster1)
df_dend$cluster1 <- as.numeric(df_dend$cluster1)
df_dend$cluster1 <- ifelse(df_dend$cluster1 == 2, 0,
                           df_dend$cluster1)

df_dend$cluster2 <- dplyr::recode_factor(df_dend$dend.clust,
                                         "Cluster 2" = "1",
                                         .default = "0")
table(df_dend$cluster2)
df_dend$cluster2 <- as.numeric(df_dend$cluster2)
df_dend$cluster2 <- ifelse(df_dend$cluster2 == 2, 0,
                           df_dend$cluster2)

df_dend$cluster3 <- dplyr::recode_factor(df_dend$dend.clust,
                                         "Cluster 3" = "1",
                                         .default = "0")
table(df_dend$cluster3)
df_dend$cluster3 <- as.numeric(df_dend$cluster3)
df_dend$cluster3 <- ifelse(df_dend$cluster3 == 2, 0,
                           df_dend$cluster3)

df_dend$cluster4 <- dplyr::recode_factor(df_dend$dend.clust,
                                         "Cluster 4" = "1",
                                         .default = "0")
table(df_dend$cluster4)
df_dend$cluster4 <- as.numeric(df_dend$cluster4)
df_dend$cluster4 <- ifelse(df_dend$cluster4 == 2, 0,
                           df_dend$cluster4)

###Create variable of groups of countries
table(df_dend$country)
df_dend$inst_design <- dplyr::recode_factor(df_dend$country,
                                            "Cyprus" = "Anglo",
                                            "Ireland" = "Anglo",
                                            "Malta" = "Anglo",
                                            "United Kingdom" = "Anglo",
                                            "Austria" = "Continental",
                                            "Belgium" = "Continental",
                                            "France" = "Continental",
                                            "Germany" = "Continental",
                                            "Luxemburg" = "Continental",
                                            "Netherlands" = "Continental",
                                            "Bulgaria" = "Eastern Europe",
                                            "Croatia" = "Eastern Europe",
                                            "Czech Republic" = "Eastern Europe",
                                            "Estonia" = "Eastern Europe",
                                            "Hungary" = "Eastern Europe",
                                            "Lithuania" = "Eastern Europe",
                                            "Latvia" = "Eastern Europe",
                                            "Poland" = "Eastern Europe",
                                            "Romania" = "Eastern Europe",
                                            "Slovenia" = "Eastern Europe",
                                            "Slovakia" = "Eastern Europe",
                                            "Denmark" = "Nordic",
                                            "Finland" = "Nordic",
                                            "Norway" = "Nordic",
                                            "Sweden" = "Nordic",
                                            "Greece" = "Southern Europe",
                                            "Italy" = "Southern Europe",
                                            "Portugal" = "Southern Europe",
                                            "Spain" = "Southern Europe") 
table(df_dend$inst_design)

df_dend$inst_design <- factor(df_dend$inst_design, 
                              levels = c("Anglo", "Continental", 
                                         "Eastern Europe", "Nordic", 
                                         "Southern Europe"))

###donwload dataframe "df_consumers.xlsx" with contry-level data
df_consumers <- read_xlsx("df_consumers.xlsx")

df_dend <- merge(df_dend, df_consumers[, c("country", "privates_respect", 
                               "trust_public", "trust_ngos","dispute_outcourts",
                               "dispute_courts","noprobs_servs","trust_justice",
                               "trust_pubadm", "trust_natgov","trust_subnatgov",
                               "trust_EU","gdp_growth","gdp_percap","trade_gdp",
                               "tax_rev","inflation","inflation_gdp",
                               "eu_membership","euro","regime")], 
                               by = "country", all.x = TRUE)

mc1 <- lm(cluster1 ~ inst_design + gdp_growth + gdp_percap + trade_gdp + 
                     tax_rev + inflation_gdp + euro + eu_membership + regime,  
                     data = df_dend)
summary(mc1)

mc2 <- lm(cluster2 ~ inst_design + gdp_growth + gdp_percap + trade_gdp + 
            tax_rev + inflation_gdp + euro + eu_membership + regime,  
          data = df_dend)
summary(mc2)

mc3 <- lm(cluster3 ~ inst_design + gdp_growth + gdp_percap + trade_gdp + 
            tax_rev + inflation_gdp + euro + eu_membership + regime,  
          data = df_dend)
summary(mc3)

mc4 <- lm(cluster4 ~ inst_design + gdp_growth + gdp_percap + trade_gdp + 
            tax_rev + inflation_gdp + euro + eu_membership + regime,  
          data = df_dend)
summary(mc4)

library("modelsummary")

models <- list("Cluster 1" = mc1,
               "Cluster 2" = mc2,
               "Cluster 3" = mc3,
               "Cluster 4" = mc4)

mdf = modelplot(models, draw = FALSE)
mdf$ID <- 1:52 
mdf <- mdf[c(5:20),]

mdf_1 <- modelplot(models, draw = FALSE, 
                   conf_level = 0.99)
mdf_1$ID <- 1:52 
mdf_1 <- mdf_1[c(5:20),]

colnames(mdf_1)[colnames(mdf_1) %in% c("conf.low", 
                                       "conf.high")] <- c("conf.low_0.99", "conf.high_0.99")

mdf_2 <- modelplot(models, draw = FALSE, 
                   conf_level = 0.90)
mdf_2$ID <- 1:52 
mdf_2 <- mdf_2[c(5:20),]

colnames(mdf_2)[colnames(mdf_2) %in% c("conf.low", 
                                       "conf.high")] <- c("conf.low_0.90", "conf.high_0.90")

mdf <- merge(mdf, mdf_1[, c("ID", "conf.low_0.99", "conf.high_0.99")], 
             by = "ID", all.x = TRUE)

mdf <- merge(mdf, mdf_2[, c("ID", "conf.low_0.90", "conf.high_0.90")], 
             by = "ID", all.x = TRUE)

table(mdf$term)
mdf$term <- dplyr::recode_factor(mdf$term, 
                                 "inst_designContinental" = "Continental",
                                 "inst_designEastern Europe" = "Eastern Europe",
                                 "inst_designNordic" = "Nordic",
                                 "inst_designSouthern Europe" = "Southern Europe")
table(mdf$term)

mdf$term <- factor(mdf$term, levels =  c("Continental",
                                         "Eastern Europe",
                                         "Nordic",
                                         "Southern Europe"))

##Figure 9

ggplot(mdf, aes(x = term, y = estimate)) +
  geom_point(position = position_dodge(width = 0.4), size = 3.5) + 
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), 
                 position = position_dodge(width = 0.4), linewidth = 1.5, lwd = 1) +
  geom_linerange(aes(ymin = conf.low_0.99, ymax = conf.high_0.99), 
                 position = position_dodge(width = 0.4), linewidth = 0.5, lwd = 1) +
  geom_linerange(aes(ymin = conf.low_0.90, ymax = conf.high_0.90), 
                 position = position_dodge(width = 0.4), linewidth = 2.5, lwd = 1) +
  geom_hline(yintercept = 0, colour = gray(1/4), lty = 2) +
  labs(x = "Predictors", y = "Coefficients") +
  theme_minimal() +
  theme(text = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.title.x = element_text(size = 15)) +
  facet_wrap(~model, ncol = 2, scales = "free_y") 

ggsave("Models1+Controls.png", width = 10, height = 8, dpi = 300)
